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ABSTRACT: Diabetic retina is a result of increase in the blood glucose level leading to damage of 
retinal blood vessels causing temporary or permanent blindness. The increase in number of cases of 
diabetic retina emphasizes the need of early detection before it reaches the severe stage. The diabetic 
retina unless other type of medical scans has multiple clots of bloods inside the vessels. Hence single 
curve is not sufficient to detect the multiple clots. A methodology is proposed here for the detection of 
multiple clots for solving a non unique solution level set segmentation. 



I. INTRODUCTION 

Segmentation [1] is defined as the portioning of the image based on certain criterions as level 
of intensity, spread of the RGB colour index, inhomogeneity in brightness etc. The objective behind 
the implementation of the of level set algorithm is based on the active contour methods abilities . The 
level set segmentation problem reduces to findings of the curves to enclose the region of interest' 31 . 
Hence the curves can be modeled directly using the control points. 




Figure 1 : Figure demonstrating the embedding of one closed contour relating to the other two closed 

contour of the same level. 

Simple formulation of the Level set algorithm 

The problem can be formulated using an active contour propagation method. Taking into account the 
segmentation boundary as a part of a surface where the contour level is zero. This is also known as zero level 
set. 

Let <P(x, t) - ±d ...(1) 

Here x is the position in the start or originating domain. T is the time and d is the distance between 
position x and the zero level set. d is positive, if x is outside zero , else d is negative. Curve of interest originates 
from — 0. which is marked. 

To evolve <P over time ,use chain rule 
•P t + "PXf. + "PyVf = . . (2) 
* t - (x t ,y t )74> = 0. .(3) 
Let (xt,yt) =n+s 
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Where n is vector normal to the front at point x and s is some arbitrary vector. Since n and s are defined 
over the entire domain of x , they are actually vector fields. 

<P t - (n-s). V<P = .. (4) 
<P t +n.V<P + s.V<P = 0. . (5) 

% - Vn. \V&\ -s.V<P = Q.. (6) 
Here Vn is some scalar. 

Values Vn and s can be viewed as two independent forces that evolve from the surface. The value of 
Vn will control how fast the surface will move in the normal direction . S dictates both direction and speed of 
evolution. The partial differential equations can then be solved when provided initial condition. 
4>(x,t = 0)..(7) 

Hence the problem can be defined as the initial value problem. 
Equation 1 can be solved by using a finite difference on a discrete spatial grid in domain of x. 

Let*(x,nAt)=*.v ■ ( 8 ) 
<|>«+i — $>V 

!J M iJ - K,M:v|-5 ; ,*:; = o..(9) 

=*"<-[V n \V *?.| - 5 ij( ..*?.]. At.. .(10) 

Curve speed and vector fields 
Vn and S can be implemented as 

A I G = (l^Jy )- . (11) 

V n = expf-oe | Vi G | ) ..(12) 
5 = /? A | W G |..(13) 

Ig is the input image blurred with a constant Gaussian filter. 0C is a constant for correcting the how fast 
Vn approaches 0, when edge encounters image and /? is a constant for controlling the strength of S. The 

parameters for the Gaussian filter are set at the user end. When the input image has small gradient the Vn 
approaches 1 and S approaches zero. However when an edge is encountered Vn approaches and s will point 
towards the edge. A larger EX will let the Vn approach at faster rate . On the other hand the o the Gaussian 

filter coefficient preserves more details of the image however the edges will be blurred. For default contraction, 
Vn is multiplied by -1. The values are calculated over central difference. 

jG _ jG 

V G *h= ^ 2 " -(14) 
I G - I G 



l" C |= -('?),,■ (16) 

N J 

V„ = expf-oc | Vi G |). .(17) 



I TJ t "j I 1 17/^ I 

ClWSDy = ' i4W 2 ^-.(18) 

(I V! G y \Xj = ^ ^..(19) 



hj = K(.\( lw G U)^CI( lv/ G | y ) i/ ..(20) 
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Entropy Condition 

The degree of randomness can be specified by considering 



ii 



(x,t) = -a(x)u x x(t) ..(21) 



* f-i x<0 

dx 




Figure 2: Diagram showing the characteristic curves with a gap, implying non unique solution. 

The figure 2 interprets that there is no unique solution for the above problem formulation . Instead the 
problem or the gap be filled by a number of curves. These curves can be used to fill the gap and multiple 
solutions can be used. However an additional constraint also called as entropy is added as shown 
u t (x,t)+u(x,t) u x (x,t)=eu xx (x,t) ..(23) 

For conversion of the equation into a hyperbolic' 41 one the viscosity term eu xx (x,t) is included which 
converts the given equation into a hyberbolic one with a unique solution. This characteristics will introduce a 
rarefaction of characteristics in the empty gap as shown in the figure below. 




Figure 3. Characteristic curve with refraction wave implying unique solution. 



Another scheme which can be implemented is the Upwind scheme. When implementing special 
attention has to be paid on how to calculate . For simplifying the problem the one dimensional problem can 

be considered 

= — ^- 



D°U. = — — . . (24) 



2 Ax 

d;U: = 

* L Ax 

D X U, = — — 

* 1 Ax 



(25) 
(26) 



These are known as central difference, forward difference and backward difference respectively. In 
order to select any one of these let us consider one dimensional problem. 
U,(x,t)- U x (x,t) =0..(27) 

U(x,0)=f(x) ..(28) 
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Where the exact solution is u(x,t)=f(x-t)=u(x-t, 0) and characteristics curves have the slope 1. Now 
considering a point A on the characteristics curve q . The value at A can be tracked back along q to a point B , 
where t=0. Hence B is in domain of dependence of Ai.e A depends on B. The information flow from A to B 
during the curve propagation. 

So an approximate scheme has to be used that propagates the information in that corresponding 
direction. In the illustration used information flows from the left to right. Therefore we use backward difference. 
In general information may flow from the left to right or from the backward to forward direction if needed. 



In order to send the information on both sides i.e forward and backward the following equations has to 



be used 

n +1 
u, - 



At 



- -(max(a^O) D x u" + min(a ; , 0) D*u"_ ± ).. (29) 



Or equivalently it can be written as 



< +1 - u" 



lad 



« + i - 2< + <_ ± ) ..(30) 



where 

it" — u(i,n)& 
a, - a(i) 




Figure 4 : Characteristic curve showing that B is in domain of A . 

This approach is known as the upwind scheme. Equation 4 is second order term that corresponds to the 
viscosity term mentioned. Using these upwind characteristics of wave propagation s and V0 can be calculated 
for every grid point. 
s.V0 = s. (0 x ,0 y ) 

= P [max(0, | V l G IJ D~ - [min(0, | V l G \ x ) D* - [max(0, | V l G | y ) D* 

- [min(0 J |7/ G | y ) D~ 0]]]].. (31) 

V n |V0| can becalculated in similar way as shown in the next step. 

V n |V0|=V0f +01-. (32) 
Where 

0! = max(O,7ji>-0- mm(0, Vj D~ .0 
: = max(O,FjD,T0 - niin(0., Vj D~ .0 



In order to preserve the stability of the system while solving the Courants- Friedrichs- Lewy (CFL 
condition) is applied. The CFL condition takes into consideration the numerical domain of dependence as shown 
in the figure below. In the particular case the analytical domain of dependence and the dashed line be numerical 
domain of dependence. 
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Figure 5. Let the solid line represent the analytical domain of dependence and the dashed line 

denote the numerical domain of dependence .In the left figure At = 1 At — 2 on the right. On the left the 

numerical domain of dependence includes the analytical domain of dependence, satisfying the CFL condition 
.However the CFL condition is violated on the right. 

In the first figure the At=l and At— 2 for the second figure. In case of figure 2 the numerical domain of 
dependence includes the analytical domain of dependence , satisfying the CFL condition. The CFL condition is 
violated in the second figure. 

In case of a CFL condition each point is updated by only its neighboringpoints information's should 
not be moving faster than one grid point per time step. To prevent the appearance of instabilities, the time step is 
implemented as 

Af = ri n tt— — .. (32) 



Where 0<c<l is some scaling constant set by the user. 
The proposed method allows the propagation of the curves in multiple directions by splitting up into two or 
more curves. The splitted curves can join or merge each other based on the propagation criteria set by the user. 

II. RESULTS 

The method of dual curve propagation was applied on various types of retinal images and the result 
was analyzed for detection of the exudate formation. The various stages of the diabetic retinopathy are analyzed. 
The multi wave propagation is used in order to search for multiple exudate formation. 



Tab 


e I. Input and Exudate Extraction and f 


stage of Diagnosis 


SI No 




Input Retinal 
Image 




Exudate 
Extraction 


Stage of Diagnosis 


1. 








Background 
Diabetic 

Retinopathy (BDR) 


2. 




R 


- 


Proliferative 
Diabetic 

Retinopathy (PDR) 


3. 






( v) 


Severe Diabetic 
Retinopathy (SDR) 



Table 1 shows three stages of diabetic retina proliferation in human eye. The stage 1 is Background 
Diabetic retinopathy(BDR) where the retinal vessels develop the initial blood clot formation. The stage 2 shows 
the second stage Proliferative diabetic Retinopathy (PDR) where the blood vessels of the retina starts developing 
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the swelling. The last stage is severe Diabetic Retinopathy (SDR) where the entire retinas and the other blood 
vessels suffers severe diabetic retinopathy and ultimately leads to blindness . 

III. CONCLUSION 

The level set segmentation of multiple wave propagation segments the retinal images and also detects 
the multiple exudate formation. The proposed method hence helps in early diagnosis of the diabetic retina , 
when the initial retinal blood clots starts to form in the retina. The proposed method is independent of any kind 
of inhomogeneity's in the input image. 
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